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ABSTRACT 

We revisit the relativistic restricted two-body problem with spin employing a pertur- 
bation scheme based on Lie series. Starting from a post-Newtonian expansion of the 
field equations, we develop a first-order secular theory that reproduces well-known 
relativistic effects such as the precession of the pericentre and the Lense-Thirring and 
geodetic effects. Additionally, our theory takes into full account the complex inter- 
play between the various relativistic effects, and provides a new explicit solution of 
the averaged equations of motion in terms of elliptic functions. Our analysis reveals 
the presence of particular configurations for which non-periodical behaviour can arise. 
The application of our results to real astrodynamical systems (such as Mercury-like 
and pulsar planets) highlights the contribution of relativistic effects to the long-term 
evolution of the spin and orbit of the secondary body. 
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The General Theory of Relativity (GR) has been one 
of the greatest achievements of the XX century. Its formu- 
lation has revolutionised our way to understand the phys- 
ical Universe and has led to the (sometimes unexpected) 
opening of entirely new research fields. One striking exam- 
ple is the formulation of a consistent theory describing the 
global behaviour of the Universe, which was impossible in 
the framework of Newtonian mechanics (Peebles 1993). 

In Einstein's theory the spacetime is thought as a four 
dimensional pseudo-Riemannian manifold whose geometry 
interacts non-linearly with matter. Because of this feature, 
in general calculations in GR require much more work than 
the corresponding ones in Newtonian gravity and this fact 
has made more difficult the understanding of the physics of 
Einsteinian gravity. 

Another consequence of this difficulty falls in the con- 
text of testing Einstein's gravity: our inability to solve in 
general the gravitational field equations limits the possibil- 
ity of devising full tests of GR. Instead, most of the in- 
vestigation on the validity of this theory was performed in 
a regime of weak field leaving the strong field regime al- 
most completely unexploredf. In fact, all the tests of GR 
performed up to now belong to the class of weak field ex- 

* E-mail: bluescarni@gmail.com 

f There are in truth other reasons why all our tests are focused 
on weak field: we live in a weak gravitational field and it is impos- 
sible to control gravitational fields like we do, for example, with 
electromagnetism. 



periments, from the so-called "classical tests" of Relativity 
(Peebles 1993) to the latest tests based, e.g., on the motion 
of interplanetary probes (Hoes et al. 2011). 

Einstein himself was well aware of the fact that the un- 
derstanding of his theory in the perturbative regime was 
crucial for its testing and its application to the problems 
in which Newtonian gravity was most successful. For this 
reason in the 1930s he developed a perturbative approach 
able to describe the subtle changes induced on the New- 
tonian evolution of the bodies in the Solar System by Gen- 
eral Relativity. The so-called Einstein-Infeld-Hoffman (EIH) 
equations are the results at first post-Newtonian level of 
this attempt (Einstein et al. 1938)J. The derivation of these 
equations was the the birth of a new field called "Relativis- 
tic Celestial Mechanics" (ROM) (Brumberg 1991; Kopeikin 
et al, 2011). 

The research in RCM has not stopped since then. The 
post-Newtonian equations have been analysed in many dif- 
ferent ways and new and more powerful approaches to RCM 
have been proposed (see Damour et al. 1991a, 1993, 1994, 
19911)). In spite of these advances, the resolution of the 
post-Newtonian equations still remains a formidable task. 
In particular, the interaction between orbital and rotational 
motion has not been fully solved analytically to this date. 

In this paper we will focus on the post-Newtonian equa- 

X To be precise, in the same period other researchers proposed 
alternative methods based on different assumptions (lock 1961; 
Papapetrou 1974), but they all lead to the same results. 
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tions for a restricted two spinning bodies system at the IPN 
level. Our target is to obtain a complete description of the 
leading correction of GR on Newtonian mechanics with a 
specific focus on the spin interaction effects. Such target will 
be accomplished using modern perturbative methods based 
on Lie series (fiori lOtiCi; Deprit 1969). This approach allows 
us to derive a fully-analytical description of the long-term 
dynamical evolution of the system using semi-automatised 
computer algebra. Advantages of the the Lie series proce- 
dure over traditional perturbation techniques include the 
availability of explicit formula to convert between the aver- 
aged and non-averaged orbital elements, and the possibility 
to straightforwardly extend the treatment to higher post- 
Newtonian orders. Using the Lie series technique, we will 
re-derive the standard precession effects typical of the re- 
stricted two-body problem in General Relativity, but we will 
also give new analytical relations describing the evolution of 
spin and orbital angular momentum of the secondary body. 
In fact, we will be able to explore the full evolution of the 
orbital parameters of the secondary body and to quantify 
the effect of the spin-orbit and spin-spin interactions. 

The use of the Lie series technique in ROM is relatively 
sparse, and, as far as we have been able to verify, it consti- 
tutes by itself a novelty in the analysis of relativistic spin- 
spin and spin-orbit interactions§. In this first paper we have 
thus chosen to focus on the restricted problem in order to 
lay the foundations of our method in the simplest case of 
interest. In subsequent publications, we will tackle the full 
post-Newtonian two-body problem with spin building on the 
formalism introduced in the present paper. 



1 HAMILTONIAN FORMULATION 

The starting point of our derivation is the well-known IPN 
Hamiltonian of the two-body problem with spin, which, af- 
ter reduction to the centre-of-mass frame, reads (Barker <k 
O'Conncll 1970, 1979; Daniour 2001) 



(1) 



Here e = is chosen as the "smallness parameter" of our 
perturbation theory and 



2 /i 2 /a 2^ r 



(2) 



is the Newtonian Hamiltonian (representing the unper- 
turbed problem). Hi expands as 



Hi = m 



-Hso+ns 



(3) 



where 



gM 



(3 + 0^ 



2r2 



(4) 



is the post-Newtonian orbital Hamiltonian, 
3 m2 \ , 3 mi 



■Hso 



2g 



1 + 



4 mi 



Ji+ 1 



4 m2 



represents the spin-orbit interaction, and 

•Hss = [3 ( Ji ■ n) ( Ja • n) - Ji • J2] 



(r X p) 

(5) 

(6) 



represents the spin-spin coupling. In these formulae, g is the 
universal gravitational constant, mi and ma are the masses 
of the two bodies, M — mi + ma, jJ. — mima/M is the 
reduced mass, u = fi/M, p = pi — — pa is the canonical 
momentum of body 1, and r = rn is the vector connecting 
body 2 to body 1. The structure of this Hamiltonian and, 
above all, the nature of the spin vectors, have been a matter 
of debate in the past few years (see Damour ct al. 2008). 
More recently, Wu & Xio (2010) have proposed a symplectic 
formulation of the spin vectors in terms of cylindrical-like 
coordinates. Here, we follow the interpretation of Barker ck' 
O'Comioll (1975, 1976, 1979) and Wcx (199.",) in identifying 
the spins Ji as the rotational angular momenta of spherical 
rigid bodies, so that 



(7) 



where li is the moment of inertia and LJi the rotational an- 
gular velocity vector of body i. 

In the present formulation of the Hamiltonian, we have 
dropped from Hss the spin-quadratic contributions due 
to monopole-quadrupole interaction'^ to focus on the lead- 
ing relativistic effects. Besides, as pointed out by Daniour 
(2001), at this PN order we can expect to obtain physically 
reliable results in case of moderate spins and orbital veloci- 
ties. 

The equations of motion generated by (1) can be ob- 
tained via the familiar Hamiltonian canonical equations. 
The orbital momenta and coordinates are represented re- 
spectively by p and r; the spin coordinates and momenta 
are represented by the Euler angles and their conjugate gen- 
eralised momenta in terms of which the components of the 
angular velocities c^i are expressed. At this stage, we are 
not concerned with the exact functional form of such depen- 
dence, as we will introduce a more convenient set of variables 
in §2. We refer to (iurfil ct al. (2007) for a derivation of the 
Hamiltonian formulation of rigid-body dynamics. 

The reduction of Hamiltonians (2)-(6) to the restricted 
case in which ma ^ mi and | Ja| ^ | Ji| reads 

Ijf p'i Smima 

2 h 2mi r 

f 1 pi 3 C/ma pi g^ml 
Wpn = mi (^---^ - - — + ^ 



„ , 25 / 3 ma , , , , , , 
Hso = -. — Ji + Ja • (r X pi) , 
r'^ V4mi ' 



(8) 



[3 ( Ji • n) ( Ja • n) - Ji ■ Ja] , 



where Ja is now considered as a constant of motion that 
can be dropped from Hn. Without loss of generality, we can 



§ We were able to locate two papers in the existing literature 
that tackle the solution of the post-Newtonian equations using 
Lie series (Kicliardsoii Kolly 19SS; Hciuibcrgcr et al. 1990), 
none of them considering spin interactions. 



^ The monopole-quadrupole interaction terms for compact bod- 
ies will depend on their internal structure and equations of state. 
For black holes, they are given, e.g., in (Damour 2001). 
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orient the reference system (now centred on body 2) in such 
a way that the constant J2 is ahgned to the positive z axis, 
so that 



J2 = (0,0, J2), 



(9) 



with J2 



2 ACTION-ANGLE VARIABLES 

In preparation for the appHcation of the Lie series peturba- 
tive approach in the study of Hamiltonian (8), we want to 
express the Hamiltonian in a coordinate system of action- 
angle variables (Arnold 1980) for the unperturbed problem. 
For the orbital variables, a common choice is the set of De- 
launay elements (AlorbidcUi 2002). For the rotational mo- 
tion, a suitable choice in our case is that of Serret-Andoyer 
(SA) variables (Gurfil ct al. 2007). Both sets of variables are 
introduced formally as canonical transformations. 

Before proceeding, we need to discuss briefly the nature 
of the canonical momenta appearing in Hamiltonian (8). 
Post-Newtonian Hamiltonians are derived from Lagrangians 
in which the generalised velocities appear not only in the ki- 
netic terms of the Newtonian portion of the Lagrangian, but 
also in the relativistic perturbation (see, e.g.. Landau & Lifs- 
chits (1975) and Straumann (1984)). As a consequence, after 
switching to the Hamiltonian formulation via the usual Leg- 
endre transformation procedure, the post-Newtonian canon- 
ical Hamiltonian momenta will differ from the Newtonian 
ones by terms of order This discrepancy will carry 

over to any subsequent canonical transformation, including 
the introduction of Delaunay and SA elements. Richardson 
& Kelly (1988) and Heimberger ct al. (1990) analyse in de- 
tail the connection between Newtonian and post-Newtonian 
Delaunay orbital elements. 

However, in the present work we are concerned with 
the secular variations of orbital and rotational elements, 
for which the discrepancy discussed above is of little con- 
sequence. For if a secular motion (e.g., a precession of the 
node) results from our analysis in terms of post-Newtonian 
elements, it will affect within an accuracy of 1/c^ also the 
Newtonian elements^ . Indeed, we will show how our analysis 
in terms of post-Newtonian elements reproduces exactly the 
formulae of known secular relativistic effects. 



2.1 Delaunay elements 

The Delaunay arguments (L, G, H, I, g, h) can be introduced 
via the following standard relations (Morbidclli 2002): 



L 

G ■- 
H ■- 



77120, 



I = M, 



G cos i, 



(10) 



n. 



Here a, e, i, M, to and Q are the classical Keplerian orbital 
elements describing the trajectory of the secondary body 
mi around the primary m^. The Keplerian elements are in 

^ Vice versa, short-term periodic variations of the post- 
Newtonian elements will have amplitudes of order and the 
precise connection with the Newtonian elements therefore be- 
comes important. 



turn related to the cartesian orbital momentum pi and posi- 
tion r via well-known relations (e.g., see .\huray k- Dorniott 
(2000)). In eqs. (10), {L,G,H) play the role of generalised 
momenta, {l,g,h) are the generalised coordinates. 



2.2 Serret-Andoyer variables 

The Serret-Andoyer (SA) variables describe the rotational 
motion of a rigid body in terms of orientation angles and 
rotational angular momentum. In our specific case, the bod- 
ies are perfectly spherical and thus the introduction of SA 
elements is simplified with respect to the general case. We 
refer to Gurfil ct al. (2007) for an exhaustive treatment of 
the SA formalism in the general case. 

We start by expressing J\ in a coordinate system where 
Ji itself is aligned to the z axis, so that Ji = (0,0, |Ji|)^ In 
order to express the components of Ji in the centre-of-mass 
reference system, we compose two rotations: the first one 
around the x axis by — / (the inclination of Ji), the second 
one around the new z axis hj —h (the nodal angle of J\, 
analogous to the h angle in the Delaunay elements set or to 
the longitude of the node in the set of Keplerian elements), 
resulting in the composite rotation matrix 



cos h — sin h cos / 
sin h cos / cos h 
sin/ 



sin / sin h 
- sin I cos h 
cos I 



(11) 



In the SA formalism, cos/ = Ji,z/\J\\ (where Ji,z is the 
z component of the rotational angular momentum in the 
centre-of-mass reference system). Thus, in the centre-of- 
mass reference system, the components of Ji read 



Ji = 



Ji| — Jj^j sin/i, — 1/1 Ji| —Jl^cosh,J\, 



(12) 

In the SA variables set, (| Ji| , Ji,z) are the generalised mo- 
menta, and ^g, are the conjugate coordinates. Note that, 
due to the spherical symmetry of the bodies, the angle g con- 
jugated to the action | Ji| does not appear in the expression 
of the Hamiltonian and I Ji| is thus conserved. 



2.3 Final form of the Hamiltonian 

We are now ready to express the Hamiltonians (8) in terms 
of Delaunay and Serret-Andoyer arguments. Before proceed- 
ing to the substitutions, in order to simplify the notation we 
rescale the Hamiltonian via an extended canonical transfor- 
mation, dividing the Hamiltonian and the generalised mo- 
menta by mi. The resulting expressions for "Hn and "Hi 
in terms of the momenta I L,G, H,G, H] and coordinates 



I, g, h,g,h\ read 



Hi 



Q ma 



2L2 ' 
1 e*m| 1 2g^ml 



(13) 



8 L* 



H 

r 



L2 



1 0^2 2 
— 3C/ m2 



^ This reference system is referred to as "invariable" in Gurfil 
ct al. (2007). 
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+ 



2J2H + 



2G2 



+ 



J2H 



f 3m2^ 
\ 2 



+ 3J2 



1 G^yH 

2 G2 



3 _ HGxyGxy \ / 7 , 

2-^2 ^ 1 COS I h - /l 



^(2/ + 25) 



J^^(l-§)cos(2/ + 25 + ft-ft 



1 G xyG xy 

4 G 



1 + ^ ) cos ( 2/ + 2(7 - ft + /i 



(14) 



For convenience, we have summarised in Table 1 the mean- 
ing of the symbolic variables appearing in Hamiltonians 
(13)-(14). We have to stress how in these expressions, r, 
Gxy, Gxy and / have to be regarded as implicit functions of 
the canonical momenta and coordinates. Specifically, 

r = r{L,G,l), (15) 



G 



G. 



Gxy (G, H) , 
xy — Gxy ^G 

f = f{L,G,l). 



(16) 
(17) 
(18) 



We are now ready to analyse the Hamiltonian using the Lie 
series method. 



3 PERTURBATION THEORY VIA LIE SERIES 

The Lie series perturbative methodology (Hori 1966; Deijr.i 
1969) aims at simplifying the Hamiltonian of the problem 
via a quasi-identity canonical transformation of coordinates 
depending on a generating function x to be determined. The 
Lie series transformation reads 



(19) 



where is the Lie derivative of n-th order with generator 
X, H' is the transformed Hamiltonian, and H is the origi- 
nal Hamiltonian in which the new momenta and coordinates 
have been formally substituted. At the first order in e the Lie 
derivative degenerates to a Poisson bracket, and the trans- 
formation becomes 



H' =HN + e {{Hn, x} + Ui) +0 (e') 



(20) 



We need then to solve the homological equation (Arnold 
1989) 



(21) 



where x ^^nd K, (the new perturbed Hamiltonian resulting 
from the transformation of coordinates) have to be deter- 
mined with the goal of obtaining some form of simplification 
in K. Since the unperturbed Hamiltonian depends only on 
the two actions L and G, we have 

|Hn,x}- TPT^- ¥^^' (22) 



where the partial derivative of x with respect to g (the coor- 
dinate conjugated to G) can be set to zero as G is a constant 
of motion. The homological equation (21) then reads 

X-/ ^^{n.-K)dl. (23) 

The technical aspects of the solution of this integral, includ- 
ing the choice of K., are detailed in Appendix A. Here we 
will limit ourselves to the following considerations: 

• the integral in eq. (23) essentially represents an averag- 
ing over the mean motion I. Consequently, the new momenta 
and coordinates generated by the Lie series transformation 
are the mean counterparts of the original momenta and co- 
ordinates; 

• thanks to the functional form of Hamiltonians (13)-(14), 
the averaging procedure removes at the same time both / and 
g from the averaged Hamiltonian; 

• the integration is performed in closed form, that is, 
without resorting to Fourier- Taylor expansions in terms 
of mean anomaly and eccentricity (iJcprit 1982; Palacian 
2002). The results are thus valid also for highly-eccentric 
orbits. 

The computations involved in the averaging procedure have 
been carried out with the Piranha computer algebra sys- 
tem (Biscani 2008). As usual when operating with Lie se- 
ries transformations, from now on we will refer to the 
mean momenta and coordinates ^L' ,G' , H' ,G' , H'^ and 

(j' , g' , h' , g' ,h'^ with their original names (^L,G, H,G, 

and (j,, g, h, g, , in order to simplify the notation. 

After having determined x from eq. (23), the averaged 
Hamiltonian generated by eq. (21) reads, in terms of mean 
elements. 



H' = "Hn + e + £1 cos (h — h 



')] 



(24) 



where £0 and £1 are functions of the mean momenta only. In 
order to further reduce the number of degrees of freedom, we 
perform a final canonical transformation^ that compresses 
the coordinates h and ft in a single coordinate ft» (replacing 
ft) and introduces a new momentum (replacing H) via 
the relations 



hf — h — h. 



(25) 
(26) 



After this transformation, the final averaged Hamiltonian 
reads 

n' ^'HN + e{To+ Ti cos ft. ) , (27) 

where J^c and J^i are functions of the mean momenta only, 

1 Jzg^Jy.ma^ ^H^j2gSn^ 15£Sn/ 
" ^ 2 G3L3 2 G5L3 8 



+ 



G^L^ 



G^'L'i 



2 G3L3 



3 HG^H,m2^ _ g'^ma^ 
2 G»L3 GL3 ' 



(28) 



dL dl 



L3 dl 



^ This transformation is reminiscent of the treatment of single- 
resonance dynamics in the N-body problem (see Lemma 4 in Mor- 
bidclU & Giorgilli (1993)). 
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Table 1. Summary and explanation of the symbolic variables appearing in Hamiltonians (13)-(11). 



Symbol 


Meaning 


Alternate expression 


L 


normalised square root of the semi-major axis 




G 


norm of the orbital angular momentum r X pi normalised by mi 


LVl - e2 


H 


z component of the orbital angular momentum normalised by mi 


G cos i 


G 


norm of spin J\ normalised by mi 


\J\\ lm\ 


H 


2 component of spin J\ normalised by mi 


G cos / 


I 


mean anomaly 


M 


9 


argument of pericentre 


UJ 


h 


longitude of the ascending node 


Q. 


h 


nodal angle of spin Ji 




f 


true anomaly 




Xi 


inverse of the moment of inertia of body 1 normalised by mi 


mi/h 


Gxy 


norm of the projection of r X pi on the xy plane normalised by mi 


^JG'^ - 


Gxy 


norm of the projection of spin Ji on the xy plane normalised by mi 


\/& - _H"2 


J2 


norm of spin J2 


\J2\ 


r 


distance between body 1 and body 2 




mi, m2 


masses of the two bodies 




g 


universal gravitational constant 





_ 3 GxyH J2Q'^Gxy*'rn2'^ 3 GxyQ^Gxy*'rn2'^ , , 



and Gxy* is the mean Gxy (see Table 1) expressed in terms 
of the new mean momentum H, : 



and 



Gxy* = \l G^ - (^H-, - hJ 



(30) 



We are now ready to proceed to an analysis of the averaged 
Hamiltonian (27). 



4 ANALYSIS OF THE AVERAGED 
HAMILTONIAN 

Preliminary considerations 

The one degree-of-freedom averaged Hamiltonian (27) is ex- 
pressed in terms of the mean momenta ^L, G, H, G, H*^ and 

conjugate mean coordinates (j., g, h*,g, . The equations of 
motion generated by Hamiltonian (27) read: 



dL 
~dt 
dG 
~dt 
dH 
dt 
dG 
~dt 
dH* 
dt 



0, 
0, 

: tJ-\ sin h* , 
0, 
0, 



(31) 



dl 
dt 
dg 
dt 
dh* 
~dt 

dg 

dt 
dh 
lit 



n2 2 



L3 



+ e 



dJ^o , 9J"i , 



■ IiG + e^^j- cos h*, 

dG 

■ e — = 1 — cos h* 

\dH, dH* 



(32) 



A first straightforward consequence of the form of the aver- 
aged Hamiltonian is the conservation of all mean momenta 
apart from H. In particular, the conservation of the mean 
momentum H* = H + H corresponds to the conservation of 
the z component of the total mean angular momentum of 
the system. We refer the reader to Appendix B for the full 
functional form of the partial derivatives of J-q and with 
respect to the mean momenta. 

We are now going to proceed to the derivation of well- 
known relativistic effect in a series of simplified cases of 
Hamiltonian (27). 



4.1 Einstein precession 

In this simplest case, all spin-related interactions are sup- 
pressed and only the classical IPN orbital effects remain. 
Formally, the spin effects are suppressed by setting J2 = 
G = and H* = H in the equations of motion (31)-(32). 
We have to point out a formal difficulty in such a direct 
substitution, arising from the fact that when the spin G is 
suppressed the angles h* and h become undefined and equa- 



tions (B7)-(B9) have a pole in Gx 



. In order to avoid 



these problems, we can write directly the equation of motion 

^ This is analogous to the coordinates in the Delaunay set be- 
coming undefined for zero inclination/eccentricity. 
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for h = h + in the general case as 



dh dh dht 
dt dt dt 

' G3L3 



SHOW' 
2 G^LS 

2 



G5G,„L3 



3 GxyJ2Q Gxy*m2 

'2~ 



G5L3 



3 Jj-g^G^^^ma" 
2 G^G^yL^ 



cos ^» 



(33) 

and verify that this equations is regular in absence of spins: 
the poles have disappeared and the the indeterminacy of /i* 
is inconsequential as the factor of cos h, becomes zero in 
absence of spins. 

The equations of motion of the mean orbital coordinates 
in absence of spins thus become 



dl 
dt 
dg 
dt 
dh 
dt 



Q 1712 



3e 



L3 

4 

Q m2 
G^L^ 



+ 6 9 



,g m2 
GL* 



15 

2 L5 



0, 



(34) 

(35) 
(36) 



whereas all mean momenta are constants of motion. We can 
recognize in eq. (35) the classical relativistic effect of the 
pericentre precession. In Keplerian orbital elements (see con- 
version formula (10)), eq. (35) becomes: 



dt 



duj 
dt 



3m|g2 
c^ai (1 — e^) 



(37) 



which, over an averaged orbit of the secondary body of pe- 
riod 



2tt 



m2g ' 



amounts to 



(38) 



(39) 



c^a (1 — e^) 

which is the classical formula for the relativistic pericentre 
precession (Einstein 1916; Straumaini 19S4). The result for 
the perturbation on the mean mean anomaly I is in agree- 
ment with the results of Richardson & Kelly (1988) and 
Heimberger et al. (1990), also obtained with a Lie series 
technique. 

4.2 Lense-Thirring effect 

The case in which only the central body is spinning corre- 
sponds to J2 / 0, G = and H* = H. As in the preceding 
case, all mean momenta are constants of motion, and the 
equations for the mean orbital coordinates become: 

dl ^ g^ml / g^ma^ HJ2G^m2'' 15 g^ma' 
dt L3 GL4 G3L4 2 



dt V G^L'i G'^LS 

dh _ Jag^ma'^ 
dt ^ ^ G'^Li ■ 



(40) 
(41) 

(42) 



Here we can clearly recognize the effect of the interaction 
of the spin of the central body with the orbit of the (non- 
rotating) secondary body. This effect is sometimes referred 
to as Lense-Thirring precession (Thirring 1918). The equa- 
tions above show that the consequence of endowing the cen- 
tral body with a spin is that of modifying Einstein's peri- 
centre precession via the additional term 

and of introducing a precession of the mean line of nodes 
with angular velocity 

^Jag^ma^ 



(44) 



Both these two formula are in agreement with the results 
of Bogorodskii (1959) and C.'ugusi & Provcrbio (1978). 

4.3 Geodetic effect 

The case in which only the secondary body is spinning 
corresponds to Ja = and G 7^ 0. This configuration is 
clearly more complicated than the previous two cases, as 
now J^i 7^ and the mean momentum H is not a constant 
of motion any more. The equations of motion for H and its 
conjugate mean coordinate read 



dH 

lit 
dh, 
~dt 



3 GxyG Gx 

2" GUF 



*ma 



sin/i* 



Qai:i + 2 G^L'i 



3 GxyQ^H^rui" 



+ 2- 



G^L^Gx 



3 Hg^Gxy.mi' 
2 G3G™L3 



(45) 
3GxyHgSn2^ 

"2 GiL^Gxy, 

cos h-t 



(46) 

A first observation is that, since Ja = 0, in this case 
there is no preferred orientation for the reference system, 
and the form of the equations of motion is rotationally in- 
variant. As an immediate consequence, the fact that the 
mean momentum _ff« = H + H is constant regardless of 
the orientation of the reference system, implies that in this 
case the total mean angular momentum vector is a constant 
of motion (whereas in the general case only its z compo- 
nent is conserved). In turn, this observation, together with 
the conservation of the norms of the mean orbital and ro- 
tational angular momenta G and G, implies that the only 
possible motion for the mean orbital and rotational angular 
momentum vectors is a simultaneous precession around the 
total mean angular momentum vector. 

It is possible to verify this result and to quantify the pre- 
cessional motion via a phase space analysis. The dynamical 
system defined by equations (45)- (46) has three equilibrium 
points. The first one, 

Gi?* 

G + G' (47) 

hi"'' = 2kTV, 

with k £ Z, corresponds to a geometrical configuration in 
which the mean orbital and rotational angular momentum 
vectors are aligned and pointing in the same direction''. Sim- 

^ The geometrical interpretation of the equilibrium points can 
be verified with the help of eqs. (25)-(26) and Table 1. 
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ilarly, the second fixed point, 



4.4 The general case 



G-G' 



(48) 



hi"^ = TT + 2kn, 



corresponds to a geometrical configuration in which the 
mean orbital and rotational angular momentum vectors are 
aligned and pointing in opposite directions. It is possible to 
check via direct substitution in eqs. (32) (with the help of 
the partial derivatives in Appendix B) that, in both these 
two fixed points, the equation of motion for h collapses to 



(49) 



thus implying that the angles h and h are constants. In other 
words, the mean spin and orbital angular momentum vec- 
tors keep their mutually (anti) parallel configuration fixed 
in space. 

The third and last fixed point. 



ry2 I ^2 f^l 



/il"' = TT + 2fc7r, 



2H, 



(50) 



corresponds to a geometrical configuration in which the pro- 
jections of the mean spin and orbital angular momentum 
vectors on the xy plane are aligned, pointing in opposite 
directions and equal in magnitude. That is, the total mean 
angular momentum vector is lying on the z axis. In this fixed 
point, the equation of motion for both h and h reads 



dt 



dh 
dt 



2G^L3 



(51) 



In other words, both the orbital and rotational mean angular 
momenta of the secondary body precess around the total 
mean angular momentum vector with an angular velocity 
given by eq. (51). Since, as remarked earlier, in absence of 
primary body spin the equations of motion are rotationally 
invariant, it is always possible to orient the reference system 
so that the total mean angular momentum vector is aligned 
to the z axis; the precessional motion described by this third 
fixed point will thus describe the evolution of mean spin and 
orbital angular momentum in the general case. 

The result in eq. (51) agrees with the relatvistic effect 
known as geodetic or de Sitter precession (dc Sitter 191()). 
In literature (see, e.g.. Barker & O'Connoll (1970, 1979) and 
Schiff (1960b, a)) one finds the following formula for the pre- 
cession of the spin of a gyroscope orbiting a large (non- 
rotating) spherical mass 7712 in a circular orbit: 



(0) 



?,Qujm2 



where u) is the average orbital velocity 

^ ^ 1 
Qm2 



(52) 



(53) 



a is the semi-major axis and n is a unit vector in the direc- 
tion of the orbital angular momentum. This formula agrees 
with eq. (51) in the small spin limit (H* — > H) and for 
circular orbits (G ^ L). 



We turn now our attention to the general case in which both 
bodies are spinning. First we will determine and characterise 
the equilibrium points of the system, then we will study the 
exact solution of the equation of motion for H. 



4-4- 1 Phase space analysis 

In the general case, Hamilton's equations for H and h, are 

dH 

dt 
dh, 
~dt 



— eJ-\ sin h-t , 



dH + dW ''''''' 



(54) 
(55) 



where J-q and J-i are given in eqs. (28)-(29) and the par- 
tial derivatives are available in Appendix B. It is clear that 
equilibrium points can be looked for in correspondence of 
either J^i = or sin ft* — 0. We will consider these two cases 
separately. 

4.4.1.1 Equilibria with J"i = The condition J^i = 
implies the existence of the fixed point 



where G'^y and G^y 



cos ft. ' = 



G^ 
J2 

-12 



H, 



(56) 
(57) 



are Gxy and Gxy* evaluated in H 
and we have defined 



J2 



m2 



(58) 



for notational convenience. The existence of this fixed point 
is subject to three constraints: 

• since > and ^2 > by definition, H'^'^^ must also 
be positive; 

• since H ^ G hy definition (as H is the z component of 
the mean specific orbital angular momentum, whose magni- 
tude is G), the equilibrium can exist only if G ^ ^2; 

• finally, since cos/il"^' must assume values in the [—1, 1] 
interval, eq. (57) implies constraints on the values of the 
mean momenta. 

In order to characterise this equilibrium point, we can com- 
pute the eigenvalues A of the Jacobian in the equilibrium 
point, and obtain 



A : 



i 2^^!^^ Gxyt,Gxy.J2 sin hi \ 



(59) 



The eigenvalues are real and of opposite signs, and the fixed 
point is thus a saddle' . 

Note that, due to the structure of eq. (54), the initial 
condition H = H^''^ results in dH/dt = regardless of the 
initial value of ft*. As a consequence, the line H = i?*-*^' is 



^ This is of course an expected result, since equilibrium points in 
one degree-of-freedom Hamiltonian systems must either be cen- 
ters, saddles or have a double zero eigenvalue (Arnold 1989). 
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an invariant submanifold of the phase space. We can char- 
acterise explicitly the time behaviour of /i* on this subman- 
ifold. For an arbitrary initial value of h,, we can solve eq. 
(55) by quadrature and obtain 



h, it) 
where 



2 tan 



A ■■ 



,^tanhf(^-^)^^^ 
A + B V 2 



3 mjC/ 



2 \Ji 



+ J2-H, 



_ 3 -(^) (^) 



(60) 



(61) 



(62) 



Note that B is non-negative by definition. We can distin- 
guish different behaviours on the line H — H^"^^ depending 
on the values of A and B. 

In the first case, B > \A\ and the radicands in eq. (60) 
are positive, so that 



lim (t) 



-2 tan 



Noting from eqs. (57), (61) and (62) that 

A 



we have 



and hence 



cos hi ' 

- A^ _ 
A + B ^ 



B' 



tan 



( lim h, (t) ) — cos hf 

\t-»oo / 



(63) 



(64) 



(65) 



(66) 



That is, when B > \A\ h» (t) tends asymptotically to the 
equilibrium and we can speak of an attractor on the invari- 
ant submanifold H = 

In the second case B < |A| and the square roots in eq. 
(60) are purely imaginary. Thanks to elementary properties 
of hyperbolic functions, we can then write 

'Va^ 



it) = 2 tan 



A + B 



B2 (It- C) - -B^ 
■ tan ' 



(67) 

That is, the evolution of /i* in time is linear with the super- 
imposition of a periodic modulation. 

In the last case, |j4| = B and eq. (60) simplifies into two 
different formula depending on the sign of A. 

If A > 0, the time evolution of ft, is given by 



ht it) — 2 tan 



C- Bt 



and ht tends asymptotically to 0. 

li A < 0, the time evolution of h, is given by 

h, (t) = 2tan"^ {-Bt-C) , 



(68) 



(69) 



and ht tends asymptotically to vr. 

The analysis above clearly shows that it is possible to 
prepare the system in such a way to induce an aperiodic 
dynamical evolution. However, the substitution of values for 
G, G, Ht and ^2 typical of planetary systems (even in the 
case of pulsar planetary systems) shows that the orbits in 
correspondence of the equilibrium point have a very small 



semi-major axis which either puts the orbit inside the main 
body or in a zone in which the approximation of weak field 
starts to fail. In this respect then the equilibrium condition is 
interesting, but for a meaningful physical interpretation an 
analysis at higher PN orders and/or within the framework 
of the full two-body problem is required. 

4.4.1.2 Equilibria with sin ft, — The case sin/i* = 0, 
similarly to the secondary body spin configuration, implies 
hi"^ = kn (with fc G Z), and leads to the following equation 
of motion for h,: 

^{H,kn) = ^e^[TZ^{H) + TZ2m, (70) 



where 



TZi {H) ^ ^2 ^5;^' -2H + 3^^ +H,+ Ja, (71) 



G2 



G2 



7^2 {H) = 



H^Gxy<!l2 Gxy*GxyJ'2 HGxy* H^Gxy^JT'. 



Gxy*G^ 



G2 



Gx 



+ ■ 



GxyG^ 



HGxy ^ H:^Gxy H H:>,GxyJ^2 



Gxy* 



Gx 



(72) 



xy* ^xy* G xy^G'^ 

and the ± depends on k. The equilibrium condition becomes 
then 



7^l {H) = ±7^2 {H) . 



(73) 



Recalling that Gxy and Gxy* are functions of H, we can 
multiply both sides of eq. (73) by Gxy*Gxy, square them 
and, after straightforward but tedious calculations, obtain 
the equilibrium condition 



/6 (H) = 0, 



(74) 



where /e (H) is a polynomial of degree six in H whose coef- 
ficients are functions of the constants of motion. The equi- 
librium points for H will be zeroes of this polynomial, but 
not all zeroes constitute valid equilibrium points: non-real 
zeroes, and zeroes that do not respect the condition \H\ ^ G 
have to be discarded. 

A noteworthy particular equilibrium point is /il'^' — 
H^'^^ = 0, with iif — and G — G. This geometrical config- 
uration corresponds to a polar orbit in which the mean spin 
vector of the secondary body is identical (in both magnitude 
and orientation) to the mean orbital angular momentum 
vector. From eqs. (31) it follows immediately that, while ft, 
remains constantly zero, ft and ft rotate with constant an- 
gular velocity 

1 mjgVa 

2' • ^'^> 

The results above indicate that the phase space has at 
most six fixed points, but our inability to solve (74) in the 
general case prevents us from giving a complete general anal- 
ysis. Thus, here we limit ourselves to show, in Figure 1, a 
sample of the phase space for a specific choice of the param- 
eters, in order to give an idea of the kind of structures that 
one is likely to encounter in this analysis. 



4-4-S Exact solution for H (t) 

We now derive an analytical solution for the time variation 
of H . The availability of a closed- form solution for H (t) al- 
lows to obtain immediately the time evolution of cos ft, via 
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Figure 1. Sample of the phase space of the system (54)-(55). The 
hne connecting the points and is an invariant submanifold. 
On this submanifold, A^o is an attractor. 



inversion of the Hamiltonian (27)*^. With H (t) and cos/i* (t) 
it is then possible in principle to integrate the equations of 
motion for the remaining coordinates. Additionally, the ana- 
lytical expression will allow us to calculate exactly the period 
of the oscillatory motion of H (t) and to make quantitative 
predictions about the behaviour of real physical systems in 
§5. 

Recalling the form of the Hamiltonian (27) and com- 
bining it with the equation for H, 



—r- = eJ^i sm/it, 
at 



we have 



dH 



U^Tl-{U' -Un -eTo) 



(76) 



(77) 



° This, of course, is not possible if singular equilibrium points are 
present or if we are dealing with the indeterminate forms arising 
when the nodal angles are undefined. 



where "H' is the Hamiltonian constant (whose value is com- 
puted by substituting the initial values of the canonical vari- 
ables in eq. (27)), and with the understanding that the plus 
sign is now to be taken when tT\ and sin /i« have the same 
sign. At this point it can be easily proven by direct calcu- 
lations that the terms in and cancel out so that the 
radicand in eq. (77) is a quartic polynomial in H . Following 
Whittakcr & Watson (1927), we can then elect 



/4 {H) = ^^Tl - {n' -Hn- eToY 

= a4, + 4asH + 6a2H^ + 4aiii'^ -I- aoH'^ 



and rewrite eq. (77) as 



dx 



= dr. 



(78) 



(79) 



where H(, is the initial value of H and to the initial time. The 
coefficients of fi (H) are reproduced in full form in Appendix 
C. The left-hand side of this equation is, apart from the sign 
ambiguity, an elliptic integral in standard form. A formula 
by Weierstrass (see Whittakcr & Watson (1927), §20.6) al- 
lows to write the upper limit of the integral on the left-hand 
side as a function of the right-hand side. Specifically, if 



(80) 



1 



{,U{t)}-^dt, 

where a is any constant, then 
X (z) — a + 



2[p{z)-lJ-{a)Y-i-J,{a)fr (a) 



- (a) 



1 

24" 



+ 7^.U{a)fr (a) 



+ \//4 (a)p' (2 



(81) 



where p{z) = p (2; 32, ffs) is a Weierstrass elliptic function 
defined in terms of the invariants 



32 ~ aoai — 4aia3 + 3a2, 
gs — 000204 + 2010203 — 



4 



00 03 



0104, 



(82) 
(83) 

and the derivatives of /4 are intended with respect to the 
polynomial variable. After the removal of the sign ambiguity 
in eq. (79), detailed in Appendix D, the time evolution of H 
is thus given by 

H (t) = Ho + 1 

2 [P (t) - j-J'^ (Ho)] ' - is /4 {Ho) fr {Ho) 



:fi (Ho) 



P{t) 



24 



fi {Ho] 



+ ^/4 {Ho) .fi" {Ho) ± ^MHo)p' {t) 



(84) 



where the ± sign is chosen based on the initial signs of sin 
and J^i, and we set the initial time io = for convenience. 

It is now useful to recall some basic properties of p {t) 
(Abramowitz Stce;un 1964). First of all, since in our case 
the invariants (?2 and are defined to be purely real and t is 
also limited to real values, p {t) and H {t) become real-valued 
periodic functions. The period of p{t), p' {t) and H {t) is 
related to the invariants g2 and ga via formula involving 
elliptic integrals and the roots ei, 62 and 63 of the cubic 
equation 



4t^ - git - 33 = 0. 



(85) 
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The sign of the modular discriminant 

A = si - 27(?|, (86) 

determines the nature of the roots ei, 62 and 63 of the cubic 
equation (85) and thus also the value of the period. 

p (t) can degenerate into a non-periodic function when 
A = 0. In particular, if 32 = ffs = 0, p (t) becomes simply 



(87) 



If <72 > and gg < 0, the case A = degenerates instead to 

3ei 



p{t) = ei + 



jsinh [(3ei)5 t] |' 



P(i) = -y + 



where ei = 62 > and 63 = — 2ei. For (?2 > and (73 > 0, 
the case A = simplifies to 

2' (89) 

{sin[(|eOn]} 

where ei > and 62 = 63 < 0. When p (t) is non-periodic, 
for t ^ 00 H (t) tends to a finite value. 

In addition to the degenerate cases of p(t), H [t) can 
degenerate also via particular values of Ho. For instance, if 
fi (Ho) = fi (Ho) — 0, it is then clear from eq. (84) that 
H (t) degenerates to 

H{t) = Ho, (90) 

which corresponds to an equilibrium point for H. It is 
straightforward, although tedious if done by hand, to check 
by substitution that the equilibrium condition found in 
eq. (56) , Ho = j Ji, leads to the condition /4 (Ho) — 
fi (Ho) = 0. Another equilibrium condition that can be 
checked via substitution is the geometric configuration in 
which both the mean orbital angular momentum and the 
secondary body's mean spin are initially aligned to the spin 
of the primary, i.e.. Ho = G and //* — G + G. This ap- 
proach allows us to sidestep difficulties in verifying this equi- 
librium condition (which results in indeterminate forms due 



to Gxy 



Gx 



and ht being undefined) within the 



dynamical systems framework. 

As a special case of the general formula (84), we can 
analyse the behaviour of H (t) when J2 = (i.e., the par- 
ticular case of geodetic precession examined in §4.3). It is 
easy to check by substitution (see formulae in Appendix C) 
that in this case the coefficients ao and ai of the quartic 
polynomial fi (H) are both zero, and thus 



52 = 3a2, 



93 = 



3 

-a2, 



(91) 
(92) 



with 



1 mig'^W 15 2mlS^ 1 G^mig^li 
^ ^2" G^L^ ^ 16^ GH? ^ T G'^L'i 



15 2^2^ 3 2G^"iie* 
¥^ cue - 8^ G616 



1 mie" 



(93) 



4 G^LS 

and the modular discriminant A is also null. As explained 
earlier and depending on the sign of gs, in this configuration 
p (t) is either a simplified periodic function of the form (89) 
or a non-periodic function of the form (88). However, it is 
easy to check by substitution of H' that 

a2 = -^^'^ [g' 2Hl + 2HA + & + 



2Gxy,oGxy,0 cos ht^O^ , 



(94) 



where the zero subscript represents initial values. It is 
straightforward to verify how the quantity in the square 
brackets is the square of the magnitude M of the total mean 
angular momentum vector, and therefore it cannot be a neg- 
ative quantity. As a first consequence, 173 is always positive 
and the behaviour of H (t) is restricted to be purely peri- 
odic. Secondly, according to eq. (89) and keeping in mind 
that ei — —02, the angular velocity of the periodic motion 
of H (t) is 



16 LSG-s 



M2 



2 L3G3 



(95) 



(the factor 2 arising from the fact that in (89) the sine is 
squared), which is in agreement with the precessional angu- 
lar velocity (51) calculated for the geodetic effect'*. 



5 APPLICATION TO PHYSICAL SYSTEMS 

We are now ready to apply the machinery we have developed 
so far to some real physical systems. We will consider three 
examples of interest: a Mercury-like planet, an exoplanet in 
close orbit around a pulsar and the Earth-orbiting Gravity 
Probe B experiment. It is clear that, because of our initial 
assumptions, we will consider highly idealised systems; we 
have to stress again that our goal here is not to produce 
accurate predictions of the dynamical evolution of realistic 
physical models, but rather to highlight the role that rela- 
tivistic effects play in the long term. 

We will find that in all cases the spin-orbit and spin- 
spin couplings induce periodic oscillations of the mean or- 
bital plane and, at the same time, of the mean spin vector. 
Such effects are typically small for the orbital plane, but, 
because of the conservation of the z component of the total 
angular momentum H^, they correspond to non-negligible 
oscillations in the orientation of the mean spin. 

We have to point out that the perturbative treatment 
outlined in the previous sections produces results in terms 
of the components of the mean angular momentum vectors 
with respect to a fixed (non-rotating) centre-of-mass refer- 
ence system (rather than in terms of obliquity and relative 
orientations). Therefore, in the geometrical interpretations 
of our results, we will always be referring to absolute (as 
opposed to relative) orientation angles. 

Given the small magnitude of the quantities involved 
in the computations, in order to avoid numerical issues we 
performed all calculations using the mpmath multiprecision 
library (.Johansson ct al. 2011). 



5.1 Mercury-like planet 

In the simplified case in which the Sun is a rigid homoge- 
neous sphere rotating around a fixed axis, and a Mercury- 
like planet is the only body orbiting it, a set of possible 
initial values for the parameters of the system are displayed 
in Table 2. The parameters correspond to a Mercury-sized 



^ In formula (51), the z component of the total mean angular 
momentum coincides with A/ in formula (95), as in eq. (51) 
the total mean angular momentum is aligned to the 2 axis. 
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Parameter 


Value (SI units) 


Lo 


2.77 X 10^5 


Go 


2.71 X 10^5 


Ho 


2.69 X IQis 


Go 


2.95 X 10*^ 


Ho 


2.93 X 10"^ 


Ji 


1.12 X lO''^ 


r\ 


6.37 X 10^ 



Table 2. Initial values (in SI units) for the parameters of a sim- 
plified Sun-Mercury two-body system. 

object orbiting at 0.47 AU from the Sun on an orbit with 
eccentricity 0.2 and inclination 7°. The plane of the ecliptic 
is identified as being perpendicular to the spin axis of the 
Sun. The planet rotates with a period of 1405 h, and the 
absolute inclination of its spin vector is close to the value of 
the orbital inclination. Its radius r\ is 2439 km. 

Figure 2 shows the time evolution of the absolute axial 
and orbital inclinations of the planet for different initial val- 
ues of the angle ft* . The period of the evolution, calculated 
from the solution for H (t) in terms of elliptic functions, is 
6.03 Ma. When /i*,o = the system is almost in equilibrium, 
as the initial axial and orbital inclinations are almost equal 
and almost parallel to the Sun's spin vector. The amplitude 
of the periodic oscillation increases together with /i,,o. The 
oscillation is much wider in axial than in orbital inclination; 



this is a consequence of the fact that for this system G 2> G: 
the conservation of //, — H + H, G and G imposes that a 
small change in the z component of the mean orbital angu- 
lar momentum, H, is proportionally a much larger change 
in the z component of the mean spin vector, H. 

The main effect that can be observed in Figure 2 is 
the geodetic precession. Indeed, in this dynamical system 
the slow rotations of both the planet and the Sun minimise 
the spin-spin interactions, and effectively relegate the spin- 
orbit effects to a precession of the planet's mean spin axis 
around the mean orbital angular momentum vector. This 
can be verified by noting that the precessional rate given by 
eq. (51) yields essentially the same period of 6.03 Ma as the 
general formula in terms of elliptic functions. 

The correlation between the oscillation amplitude and 
/i,,o has a simple geometrical interpretation: when ht,o = 0, 
the mean spin and orbital angular momentum vectors share 
the same inclination (7°) and nodal angle, and therefore 
they are parallel (i.e., their relative angular separation is 
zero) and no precession motion takes place; as the difference 
in initial nodal angles (i.e., /i*,o) increases, the relative an- 
gular separation between the two vectors increases too and 
the mean spin precesses around the mean orbital angular 
momentum vector. When /i«,o = tt, the initial angular sepa- 
ration reaches the maximum possible value, and the projec- 
tion of the precession motion on the z axis (which is what 
is visualised in Figure 2) reaches its maximum oscillatory 
amplitude too. 
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Figure 3. Time evolution of the absolute axial (top) and orbital (bottom) inclinations of a pulsar planetary two-body system. In the 
bottom panel the quantity on the y axis is the difference from the initial value of orbital inclination. Time is measured in years. 



Parameter 


Value (SI units) 


Lo 


3.33 X lO^-* 


Go 


3.32 X lO^-* 


Ho 


3.12 X 10^"' 


Go 


5.32 X 10^" 


Ho 


5.23 X 10^0 


J2 


4.83 X lO"*! 


ri 


2.76 X 10^ 



Table 3. Initial values (in SI units) for the parameters of a pulsar- 
Jovian planet two-body system. 



ial and orbital inclinations of the planet for different initial 
values of the angle ft*. The period of the evolution is much 
shorter than in the previous case, amounting to 40 years 
circa. This system is lacking an almost-equilibrium point 
for ht.Q — 0, as the initial values of the inclinations are not 
close (although the amplitude of the oscillation is still cor- 
related to ft*,o). As in the previous case, the amplitude of 
the oscillation for the axial inclination dominates over the 
oscillation in orbital inclination, but since the G/G ratio is 
now larger the oscillation in orbital inclination is also larger, 
reaching almost 0.01° when htfi ~ 180°. 



5.2 Pulsar planet 

In this second case, the central body is a millisecond pulsar 
with a Jupiter-like planet in close orbit. The parameters of 
the system are displayed in Table 3, and they are similar 
to the estimated parameters of the PSR J1719-1438 system 
(Bailes ot al. 2011): the mass of the star is 1.4 M©, its spin 
period is 5.8 ms and its diameter is 20 km, while the planet 
has a mass roughly equal to Jupiter (1.02 M^^) and a radius 
of 0.4 r^^, orbiting on a moderately-inclined (i — 20°) near- 
circular (e = 0.06) orbit with semi-major axis of 600000 km. 
Lacking estimates on the rotational state of the planet, we 
hypothesize a spin period similar to Jupiter (10 h) and an 
absolute inclination of the spin vector of 10° . As in the pre- 
ceding case, the plane of the ecliptic is identified as being 
perpendicular to the spin axis of the star. 

Figure 3 shows the time evolution of the absolute ax- 



5.3 Earth-orbiting gyroscope 

In this last example, we consider the gravitomagnetic effects 
on a gyroscope in low-orbit on board of a spacecraft around 
the Earth. The parameters of the system are taken from the 
experimental setup of the Gravity Probe B mission (Evcritt 
ct al. 2011): the orbit is polar {i = 90.007°) with a semi- 
major axis of 7027 km and low eccentricity (e — 0.0014). 
The gyroscope consists of a rapidly rotating quartz sphere 
(38 mm diameter) whose spin axis is lying on the Earth's 
equatorial plane (i.e., the spin plane is also "polar"). The 
spin vector of the gyroscope and the orbital angular mo- 
mentum vector of the spacecraft, both lying on the equato- 
rial plane, are perpendicular to each other. Translated into 
mean Delaunay and Serret-Andoyer parameters, this initial 
geometric configuration implies H ~ 0, H, ~ 0, Gxy ~ G, 
Gxy* ~ G and h, = ±7r/2 (with the sign depending on 
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the values of h and h - note that in Evcritt ct al. (2011) 
the configuration shown in Fig. 1 implies = —tt/2). The 
substitution of these values in the general formula for H (t) 
yields a period of roughly 195 ka. 

Turning now our attention to the equations of motion, 
we first note that, given the timescale of the time-evolution 
for H, we can consider H and fe* as constants for the du- 
ration of the Gravity Probe B experiment. With this as- 
sumption, the equation of motion for h radically simplifies 
to 



dh _ 1 mlg^Ji 



(96) 



and thus h evolves linearly with time. This is the expected 
frame-dragging precession, which, after substitution of the 
appropriate numerical values, amounts to circa 0.04" per 
year (Ei-cnitt al. 2011). The effect corresponds to a move- 
ment of the mean spin axis of the gyroscope along the equa- 
torial plane in the direction of the rotation of the Earth. 
The other expected relativistic motion is the geodetic ef- 
fect, which is a drift of the mean spin axis on the orbital 
plane in the direction of the orbital motion. Translated into 
mean Serret-Andoyer variables, this effect affects the ele- 
ment H = Hf—H. Recalling that H, is a constant of motion 
and that in our experimental setup H ~ and Gxy ~ G, 
the time evolution of H is then given by 



dH 
~dt 



dH 
'~dt 



3 G^y.maC/* . , , 

r L3G2 (97) 



Recalling now that H = GcosI by definition, where I is 
the inclination of the mean spin vector with respect to the 
z axis and G is a constant of motion, we can write 



dH ^ d cos I 

dt ^ dt ' 



(98) 



and since G ~ Gxy* and / ~ 7!"/2, a Taylor expansion of 
cos / gives the simplified formula 



dl 3 mig-^ 
dt^rUG^''''^* 



(99) 



Depending then on the orientation of the angular momen- 
tum and spin vectors, /i, = ±7r/2 and the time variation of 
I becomes then 



dJ_ 
di 



^3 mte* 



(100) 



2 L3G2' 

Substitution of the numerical values in this formula yields 
the expected value of circa ±6.6" per year. 



6 CONCLUSIONS AND FUTURE WORK 

In this paper we have analysed the post-Newtonian evolution 
of the restricted relativistic two-body problem with spin. 
Our analysis has been performed via a modern perturbation 
scheme based on Lie series. The first important advantage of 
this method is that it allows us to work out all the calcula- 
tions in a closed analytic form, thus permitting an analysis of 
long-term perturbative effects which are otherwise very dif- 
ficult to detect with numerical techniques. Secondly, the Lie 
series formalism can be iterated to higher post-Newtonian 
orders using essentially the same procedure employed in this 



study, and the connection between the transformed coordi- 
nates at each perturbative step is given by explicit formula 
(contrary to the classical von Zeipel method (von Zeipol 
1916)). Finally, the Lie series methodology can be easily im- 
plemented on computer algebra systems, thus relegating the 
most cumbersome parts of the development of the theory to 
automatised algebraic manipulation. 

Our analysis reproduces well-known classical results 
such as the relativistic pericentre precession, and the Lens- 
Thirring and geodetic effects. In addition, we are able to 
thoroughly investigate the complex interplay between spins 
and orbit, and we provide a novel solution for the full av- 
eraged equations of motion in terms of Weierstrass elliptic 
functions. This result is particularly interesting as it estab- 
lishes a connection with recent developments in the exact 
solution of geodesic equations (see, e.g., Hackmanu ot al. 
2010; Scharf 2011; Gibbous & Vyska 2012), which also em- 
ploy elliptic and hyperelliptic functions. 

From the mathematical point of view, our investigation 
highlights the existence of a rich dynamical environment, 
with multiple sets of fixed points and aperiodic solutions 
arising with particular choices of initial conditions. From a 
physical point of view, we have shown how some of these ex- 
otic configurations appear in correspondence of initial condi- 
tions for which the approximation of our model starts to be 
unrealistic. In this sense, the generalisation of these results 
to a full two-body model (to be tackled in an upcoming 
publication), will provide a more realistic insight into the 
physics of these particular solutions. 

The application of our results to real physical sys- 
tems shows how relativistic effects can accumulate over 
time to induce substantial changes in the dynamics. In 
particular, the absolute orientation of the spin vector of 
planetary-sized bodies in the Solar System undergoes sig- 
nificant changes, driven chiefly by the geodetic effect, over 
geological timescales. In close pulsar planets, the evolution 
time shortens drastically to few years, and the effects on 
the planet's orbit increase by multiple orders of magnitude. 
In this sense, our results could be used to devise new tests 
of General Relativity based on precise measurements of the 
orbital parameters of pulsar planets and similar dynamical 
systems. 
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APPENDIX A: CLOSED-FORM AVERAGING 

We detail in this appendix the solution of eq. (23), here 
reproduced for convenience: 

X-/ ^jn.-lC)dL (Al) 

From eq. (14), the structure of Hi is as follows: 

-Hl=A+^ + ^ + ^ [A3a + AsbCOS {h - ft)] 

+ ^ \Bo cos (2/ + 2g) + Bi cos (2f + 2g + h-hJ 

+B2 cos (2f + 2g - h + hj^ , ( A2) 

where the coefficients Ai and Bi are functions of the mo- 
menta only. As previously remarked, in this expression we 
have to consider r and / as implicit functions of L, G and 
I. For the solution of eq. (Al), it is convenient to look for a 
IC in the form 



IC = Ao + ICi , 



(A3) 



where is a function of the momenta only from eq. (A2). 
Eq. (Al) thus becomes: 



X 



{Hi - Ao) dl - / K-idl 



(A4) 



We use now the differential relation (Murray & Dermott 
2000) 

2 2/^2 2 

to change the integration variable in the first integral of eq. 
(A4) from / to /: 



X = 



G 



{Hi-Ao)df 



ICidl, 



(A6) 



where the first integrand in eq. (A6) has now become 

{Hi - Ao) = rAi + A-z + ^ ^Aaa + Asb cos (h - ftj] 

+ - \Bo cos (2/ + 2g) 
r L 

+ Bi cos [2f + 2g + h- 

+B2Cos(2f + 2g-h + h)j^, (A7) 
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and the coefficients Ai and Bi are still functions of the mo- 
menta only. We now perform a further split of the first inte- 
grand in eq. (A6), separating the term in its own integral: 



X= I rAidf+^ I \a2 + ^ ^A:ia + Asb cos { h - h 



G J ' G _ 
+ ^ \Bo cos (2/ + 2g) + Bi cos (^2/ + 254- 

+ B2 cos (2f + 2g-h + h)\\df - 



h-h 



tCidl. 



(A8) 



defined as a 27r-periodic function in the angular coordinates. 
By plugging eq. (A13) into eq. (A.3), we can now see how 



K, = £0 + £1 cos Ih — h 



(A15) 



where the £ coefficients are functions of the momenta only. 
Finally, via this definition of IC and eqs. (21) and (20), we 
obtain eq. (24): 



n' ='HN + e [£q + £1 cos (h - /i) j 



(A16) 



The next step is to change the integration variable in the 
first integrand of eq. (A8) from the true anomaly / to the 
eccentric anomaly E via the standard differential relation 



df =- \J\ - e'^dE ■ 



LG 
rQm2 



dE, 



(A9) 



where the eccentricity e is to be considered an implicit func- 
tion of L and G and the eccentric anomaly E an implicit 
function of L, G and I. Eq. (A8) now reads: 



X = 



1 f AiLG 



dE + ^ j + ^ [A'ia + Asb cos 



G J Qm2 
+ ^ \Bo cos (2/ + 2g) + Bi cos (2f + 2g + 



h-h 



h-h 



+ B2Cos(2f + 2g-h + h) \ \df-^^ 



Kidl. 



(AlO) 



By substituting the standard relation 

1 1 + e cos / Qm2 , Qm2 , / a 1 1 ^ 

r = ^(TT^ = + G^ecos/, (All) 

in the second integrand of eq. (AlO), and after straightfor- 
ward algebraic passages, eq. (AlO) becomes 



APPENDIX B: PARTIAL DERIVATIVES OF To 

AND J"i 

The full form of the partial derivatives of To with respect 
to the mean momenta appearing in eqs. (32) is: 



dTo 
dL 



dTo 
dG 



dTo 
OH 



dTo 



2 G^L* 



15 £W 
2 L5 



9 Hg^H^mi'^ 
2 CPIA 



+ 

+ 9 



2 G3L4 2 G5L4 
g^ma^ _ 9 HJ2g*m2^ 

2 G^L^^ ' 

3 n rr2n4,^_4, 



GL* 
9 HJaS'^ma 
2 



G5L4 
(BI) 



G4L3 

g*m2^ _ 3 J2g*H*m2^ 
2 



QH^^gSn^ 15 H'J2g^H^m2-' 
2 G'^L^ ^ y GHJ 



15 H\J2g^m2-' 



G2L3 
9 Hg^H,m2^ 
2 



G4L3 



G6L3 



G4L3 
3 Ja^W^ 
2 G3L3 



HJ2g^H,m2-' 9H\j2gSn2 
G^L3 2 G^L3 



(B2) 

3 



_ Hg*m2* sg^jKm^ 

G^L'^ 2 G^L» ' 

1 Ja^W^ 3m^yn2^ 

2 G^L^ 2 G3L3 



■iH^J2g*m2-' 
2 



G5L3 



(B3) 
(B4) 



The partial derivatives of T\ are instead: 



+ y ^ 'Djknm COS ( jf + kg + nh + mh ) df 



02 2 I ^idl, 
g^mi 



(A12) 



where the C and V coefficients are functions of the momenta 
only, and the true anomaly / always appears in the cosines 
associated to the D coefficients. We can now choose ICi as 



GL3 



[Ci 4-C2,o +C2,icos (^ft- /i)] , (A13) 



so that the final expression for the generator x becomes 



X=-^lCi{E-l)+ [C2,o+C2,icos (ft-/^)] if -I) 



+ 1 ^ 'Djknm cos (^jf + kg + nh + mhj df \ , 

j^O,k,n,m I 

(A14) 

where the integral in this expression is trivially calculated 
(as the D coefficients do not depend upon /). x has thus been 



dTi 
dL 



9 Gxyg^Gxy*m2'^ 9 GxyHJ2g^Gx 



,m2 



G^L^ 



(B5) 



dTi _ 15 GxyHJ2g''G 
~dG ^ ~2 



2 GSL3 

, 3 g'^Gxv*'rn2'^ 
+ ^ 



2 G^GxyL^ 



9 Gxyg Gxy*m2 

2 G4L3 
3 HJ2g^Gxy,m2^ 
2 G-^GxyL^ ' 



dTl 3GxyH^J2g^m2'' , 2,H^J2g*Gxy*m2 



OH 



+ 



G'-L^Gxy* 

3 Gxyg^H*m2'' 
2 G^L^Gxv, 



+ 



G^GxyL^ 

3 GxyJ2g^Gxytrn2^ 
2 G5L3 



3 Gxylig 1712 
2 G^L^Gxy* 

ZGxyHJ2g^H,m2 
2 



3 Hg^G: 
2 

3 



*m2 

G^GxyL^ 



dTl 

dG 
dTl 
dH, 



G'-L^Gx 
3 Gxyg*^Gm2'^ 
2 G^L^Gxy, ~ ''■ 
3 Gxyg^H^m2^ 
2 G^L^Gx., 



3 GxyHJ2g^Gm2^ 
2 G^L^Gxy, ' 

^GxyH^jhg^^ 

2 G5L3G,^, 



(B6) 



(B7) 
(B8) 
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^ 3 G^yHg*m2* 
2 G^L^G^y, 



+ 2 



3 G^yHJ2Q^H,m2-' 



G^^L^G^ 



(B9) 



APPENDIX D: REMOVAL OF THE SIGN 
AMBIGUITY IN THE EXPLICIT SOLUTION 
FOR H 



APPENDIX C: COEFFICIENTS OF THE 
QUARTIC POLYNOMIAL 

Here follows the full form of the coefficients of the polyno- 
mial /4 (H) in eq. (78): 

°o 4 gio^e 4 G8L6 ' ('^^> 

_ 9 Ja^^CVma^ 15 J2Vi?*eW^ _ 45 Ja^Vm/ 
^ 8 + 8 G»L6 32 G^L^ 

9 Jag^e'ma^ 3 Jag^^iG'ema'^ , 3 Jag^ema^ 
2 C^LB 8 " 



+ 



8 G^LS 



3 Jag^-H'ema" 

4 G5L3 



(C2) 



4 G3L3 8 G6L6 8 G^L^ 



+ 



15 e^^e^ma* 1 G^U'erm* 1 g'^ema 



6 



16 G^L^ 



2 G3L3 



4 G^LS 



15J2^^ff*An/ 3 Ja^e^^G^e^ma*' 3 S^^G^e^ma* 



16 



G^L^ 



G»L^ 



G^LB 



ISOVma^ 1 Jag'^XiG^iJ^gma^ 
¥ G4L6 ^ 4 G5L3 

1 J2Q*n'H,em2^ _ 1 Jag'^iJ^ema^ 

2 G^L3 4 G5LJ5 
7 Jag^H,£^ma^ 

2 G6L6 



(C3) 



_ 3 Jag'^JiG^ema^ 3 JiG'^U'emi^ _ 45 G^EU^rn^ 
^ 8 G3L3 4 G3L3 32 G^U 

3 G*n'H.em2^ 3 JiG^H^e^np/ 

4 G3L3 ^ 4 GSLB 

3 JaOWa^ _ 3 JaVH^eW^ 3£^ff^em/ 
8 G3L5 8 G6L6 8 G^L^ 



+ 



9 JaS^e^ma^ 9 JaS^G^e^ma^ 27 S^^.e'ma^ 



+ 



04 



4 G^LS 8 G8L6 ' 8 
45 Ja^Vma^ _ 3 G^X^&H^em2^ 
32 GHJ 8 GH7^ ' 

1 J2G^IiG^H,em2^ IJiG^H^em^ 

2 G^L^ 2 G3L5 



G*L^ 



(C4) 



15 e'ema'' , ^a , 1 ^'XiG^ma' 

Y^J^^^^'^ +2 

1 G*m2^ _ g'H'ma^ _ 15 g^XiG'ema'* 
4 



6 



+ 



L2 



8 



, J2G''n'mem2-' _,,a , 15 g-^^'ema^ 1^2/^4 



1 JaV^ffJeW 
4 G6L6 



9 



G2L6 



+ 3 



GL3 



225 a^^e-'ma*' 45 e'^e'^ma*' 9 S'^G'^e^ma** 
li4 L» T GL'^ ~^ 4 G*L6 

GL5 4 G4L6 G4L6 
15 Jag^g.e'ma^ g^-H'ema^ 
8 G3L7 ^ ■ 



(C5) 



In this appendix we detail the removal of the sign ambiguity 
in eq. (79), here reproduced for convenience: 



dx 

Ho xZ/tM 



dr. 



(Dl) 



We first note how the sign change happens in correspondence 
of the zeroes of the polynomial /4 (x) . Secondly, we note 
how Weierstrass's elliptic function p (z) is an even function, 
whereas p' (z) is odd. Finally, from eq. (81), we introduce 
the shorthand 

W{a;z) ^ 



2[p(z)-^J-(a)Y^^f,{a)ft (a) 



- ^./r (a) 



+ ^/4 (a) /r (a.) + VU (a)p' (z) 



(D2) 



in order to simplify the notation. Without loss of generality, 
we suppose that the initial values of /i* and H are such 
that the initial value of the time derivative of H is positive 
and the plus sign is then chosen in the integrand. We can 
integrate from Ho to the next sign change in correspondence 
of a zero Hi — H (ti) of /4 (H) and obtain 



Hi=Ho + W{Ho;ti-to). 



(D3) 



Next, taking Hi as new initial value, and H2 = H {12) as 
the next zero of /j [H), we can write 



H2 



Hi+WiHu±tiTt2), 



(D4) 



where the plus-minus signs take into account the possibility 
of sign change of the integrand. However, since Hi is now a 
zero of /4 (-ff), the W function simplifies to 



W{Hi;ti~t2) = 



4p(ti-ta)- 



(D5) 



thus becoming an even function and allowing us to choose 
either sign freely: 



H2=Hi+W{Hv,t2-ti). 



(D6) 



We can repeat this process as many times as needed, and 
write for a generic value H 

n-l 

H = Ho + W {Ho;ti - to) + Y,^ iH,;U+i ~ ti) 

+ W{Hr,;t-~t„), (D7) 

where Hi, . . . ,H„ are zeroes of (H). This expression is 
the same we would have obtained if we had just fixed in 
eq. (Dl) the sign of the integrand once and for all based 
on the initial signs of sin/i, and J^i and repeated the same 
procedure. Analogously, when the initial sign is negative we 
can write 

H = Ho + W (Ho; to - ti) + ^ W (H,;t^ - t,+i) 

i = l 

+ W{H„;t„-t). (D8) 
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It follows then that the sign in eq. (Dl) can be selected 
according to the initial signs of sin h, and J-i , and the time 
evolution of H is thus given by eq. (84). 

This paper has been typeset from a TfrjX/ M^jX file prepared 
by the author. 
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